/*
** This program is a modified version of
** 
** Bias in Conditional and Unconditional Maximum Likelihood Estimation
**  Ethan Katz (ekatz@fas.harvard.edu)
**  last modified 02/06/01
**
**  The goal of this program is to compare the conditional and unconditional
**  maximum-likelihood estimators of fixed-effects logit models.  It generates
**  a panel data set with "n" units, "t" time periods, one structural
**  parameter set to "xcoeff", a constant set to "concoeff", and "n"-1
**  incidental parameters set to 1.  It then estimates the parameters using
**  both conditional and unconditional maximum likelihood.  It repeats this
**  process "sims" times, storing the structural parameter estimates each
**  time, and then calculates the average estimates.  These averages
**  can be taken as estimates of the expectations of the parameter estimates,
**  and their deviations from the true parameter can be taken as estimates of
**  bias.  Finally, the program calculates the Monte Carlo standard errors
**  (the error from using a finite number of samples to estimate the mean of
**  the parameter estimate sampling distribution).
**
** date: 02 02 2005
** this modification solves the error in Katz's original program as it does allow to include fixed effects
** Tom Coup?  tcoupe@eerc.kiev.ua
*/

version 8.0
set memory 125000
clear
matrix define results=J(15,16,0)
matrix define results200=J(15,16,0)
foreach time in 2 4 8 12 16{ 
foreach number in  100 { /*200 500*/
quietly {
drop _all



noisily {

local n=`number'		/* number of units */
local t=`time'			/* number of time periods */
local concoef=-1		/* coefficient on constant term */
local xcoeff=0.5	/* coefficient on X */
local sims=500		/* number of simulations */

}

set matsize 800
local nt=`n'*`t'		
local n1=`n'+1
local n2=`n'+2
local t1=`t'+1


/* this is to generate x's */ 

set obs `n'
matrix f=matuniform(`n',`t'+1)
svmat f


/* for generating x's and fixed effects you have to choose the right combination below!*/ 



forvalues i=2/`t1' {
local ii=`i'-1

generate fix`ii'=(invnorm(f1)+invnorm(f`i'))*(sqrt(2)/2)  /* this is to have real fixed effects specification with correlation between x and fixed effects */ 

*generate fix`ii'=invnorm(f`i')  /* this is to have  specification without correlation between x and fixed effects */ 

*generate fix`ii'=f`i'  /* this is to have Katz specification with uniform x without correlation between x and fixed effects */ 
}


/* this is to generate fixed effects */ 

generate fef=invnorm(f1) /* this is to have fixed effects specification with correlation between x and fixed effects */ 

*generate fef=0 /* this is to have fixed effects specification without correlation between x and fixed effects */ 

*generate fef=1 /* this is to have katz's specification (in his paper) without correlation between x and fixed effects */ 

mkmat fef,matrix(fef1) 
matrix define ff=J(1,`t',1)
matrix fefe=ff#fef1
svmat fefe

generate nn=_n

reshape long fix fefe, i(nn)   /* fix is the x values , fefe the fixed effect - we reshape it so as to get panel setup rather than cross section setup*/ 

rename fix x

gen p=	1/(1+exp(-(`concoef'+fefe+`xcoeff'*x)))



matrix define d=J(1,`n',0) /* creates matrix of unit-specific dummy variables */
svmat d
replace d1=1
set obs `nt'
local i=1
while `i'<=`n' {
	replace d`i'=0
	local j=(`i'-1)*`t'+1
	local k=`i'*`t'
	replace d`i'=1 in `j'/`k'
	local i=`i'+1
	}

/* note that we make a dummy for each individual*/


gen group=group(`n')

if `sims'>`nt' {
	set obs `sims'
	}
gen clest=.		/* stores clogit conditional estimates of X coefficient */
gen lest=.			/* stores unconditional with dummies estimates of X coefficient */
gen nfest=.			/* NEW stores unconditional without dummies (simple logit) estimates of X coefficient */

local i=1
while `i'<=`sims' {	/* draws Y from a Bernoulli distribution and estimates parameters */
	gen temp=uniform()
	gen y=0
	replace y=1 if temp<p
	clogit y d1 x, group(group)
	replace clest=_b[x] in `i'
	logit y d* x, nocon
	replace lest=_b[x] in `i'
	logit y x             /*new*/
	replace nfest=_b[x] in `i' /*new*/
	if `i'<`sims' {
		drop temp y
		}
	local i=`i'+1
	}

}


if `n'==100 { 
	if `time'==2 {
	local counter=1
	}
	else if `time'==4 {
	local counter=2
	}
	else if `time'==8 {
	local counter=3
	}
	else if `time'==12 {
	local counter=4
	}
	else if `time'==16 {
	local counter=5
	}
}
if `n'==200 { 
	if `time'==2 {
	local counter=6
	}
	else if `time'==4 {
	local counter=7
	}
	else if `time'==8 {
	local counter=8
	}
	else if `time'==12 {
	local counter=9
	}
	else if `time'==16 {
	local counter=10
	}
}

if `n'==500 { 
	if `time'==2 {
	local counter=11
	}
	else if `time'==4 {
	local counter=12
	}
	else if `time'==8 {
	local counter=13
	}
	else if `time'==12 {
	local counter=14
	}
	else if `time'==16 {
	local counter=15
	}
}

matrix results[`counter',1]=`n'
matrix results[`counter',2]=`time'

summarize clest					/* displays mean conditional estimate of X coefficient */
matrix results[`counter',3]=r(mean)
local cstd=sqrt(r(Var))/sqrt(`sims')	/* Monte Carlo standard error */
display `cstd'
matrix results[`counter',4]=`cstd'
summarize lest					/* displays mean unconditional with dummies estimate of X coefficient  */
matrix results[`counter',5]=r(mean)
local std=sqrt(r(Var))/sqrt(`sims')	/* Monte Carlo standard error */
display `std'
matrix results[`counter',6]=`std'


summarize nfest					/* NEW displays mean unconditional estimate without dummies (simple logit) of X coefficient */
matrix results[`counter',7]=r(mean)
local cstd=sqrt(r(Var))/sqrt(`sims')	/* Monte Carlo standard error */
display `cstd'
matrix results[`counter',8]=`cstd'

summarize p					/* shows how rare the event is */
matrix results[`counter',11]=r(mean)
matrix results[`counter',12]=`concoef'
matrix results[`counter',13]=`xcoeff'


matrix list results



}
}
